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Abstract 

The Oslo sandpile model, or if one wants to be precise, ricepile model, is a 
cellular automaton designed to model experiments on granular piles displaying 
self-organized criticality. We present an analytic treatment that allows the 
calculation of the transition probabilities between the different configurations 
of the system; from here, using the theory of Markov chains, we can obtain 
the stationary occupation distribution, which tell us that the phase space is 
occupied with probabilities that vary in many orders of magnitude from one 
state to another. Our results show how the complexity of this simple model is 
built as the number of elements increases, and allows, for a given system size, 
the exact calculation of the avalanche size distribution and other properties 
related to the profile of the pile. 
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I. INTRODUCTION AND DEFINITION OF THE OSLO MODEL 



Self-organized criticality (SOC), born from the deep insights of Bak et al, deals with 
the emergence of scale invariance in slowly-driven nonequilibrium systems [1,2]. The phe- 
nomenon is illustrated with the archetypical example of a pile of sand, and realized in 
computer simulations of diverse sandpile models, which are mainly based on the original 
Bak- Tang- Wiesenfeld (BTW) model [3]. However, the relevance of SOC for real granular 
matter was unclear until Frette et al. [4] performed experiments on a 1 + 1— dimensional pile 
of rice; these experiments and some others [5,6] were modeled with a cellular automaton 
introduced in Ref. [5], later called the Oslo model. 

The Oslo model has the interest of being (as far as we know) the first SOC sandpile 
model or, more appropriately, ricepile model, able to reproduce experimental results. For 
the avalanche properties [4], the concordance with experiments is only qualitative, whereas 
for the transport of individual grains [5], and for the surface roughness [6], the agreement 
is also quantitative. Moreover, the Oslo model is remarkable as a simple model of SOC, 
because it displays this nontrivial behavior in one dimension. 

The model is designed to mimic the experimental situation in Refs. [4-6]: grains are 
slowly added at a fixed position on a quasi one-dimensional substrate which is in between 
two parallel vertical plates; just at the (let us say) left of the position of addition a wall 
prevents the falling of the grains; on the other side, the right boundary is open. The model 
assumes a discrete space, x — 1,2 ... L, from left to right, as well as discrete time and field 
(the height of the pile, or number of grains). The grains pile up in columns until the local 
slope somewhere is too large, then the upper grain becomes unstable and is transferred to the 
next column to the right (from x to x + 1). This transference can induce further instabilities 
and therefore a chain reaction or avalanche. A slow-driving condition imposes that during 
the avalanches the external addition of grains is interrupted; this implies that the evolution 
of the model (and of the experiments) takes place in two separated time scales: a slow time 
scale for the grain addition and a fast time scale for the evolution of the avalanches. 
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In terms of the height of the x column, h(x), and the local slope, defined as z(x) = 
h(x) — h(x + 1), taking h(L + 1) = 0, these prescriptions are expressed in the following rules 
[5]: 



if z(x) < z th {x) \fx h(l) -> h(l) + 1, 

h(x) — > h(x) — 1, 

if z(x) > z t h{x) =>• < 



(1) 



1) -> /i(a; + 1) + 1, 
z t h(x) — > rand 

(where the update is supposed to take place in parallel). Here z t h(x) refers to a local 
threshold, which rather than being constant changes with every toppling at a; to a random 
value rand, chosen as 



rand = < 



(2) 



1 with probability p, 

2 with probability q = 1 — p. 

These simple rules, Eqs. (1) and (2), together with the boundary condition, h(L+l) = 0, 
completely define the Oslo model. Nevertheless, it is convenient to express rule (1) in terms 
only of the slope, turning out to be, 

if z(x) < z th {x) \/x =>• z(l) -»• -z(l) + 1, 

z(x) — * z(x) — 2, 



if z(x) > z t h{x) 



z(x — 1) — > z(x — 1) + 1, 
z{x + 1) -> z(x + 1) + 1, 

^th(^) — >■ rand, (3) 
for x = 1, 2 ... L — 1 and 

z(L) ^ z(L) - 1, 
z(L-l) ^ z(L-l) + l, 
Zfh(L) — > rand 

(taking into account that at rr = the z variable is not defined). Notice how the dynamics 
of the grains has allowed to define a different dynamics for the units of slope, which can be 



if z(L) > z th (L) < 
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considered as some kind of virtual particles. Both dynamics are conservative, except at the 
boundaries, which curiously are reversed: the open end for the grains at x = L is a closed 
boundary for the slopes and vice versa at x — 1. It is important to have this in mind to 
avoid confusions. 

The Oslo model is essentially the one-dimensional BTW model, but with fixed addition 
at x = 1 and an open boundary condition for the grains at x — L. The key different 
ingredient, which makes the model critical in one dimension, is the selection of dynamically 
changing thresholds, to account for the heterogeneities of a real system. In this way, whereas 
the randomness in the BTW sandpile is external, in the Oslo model it is internal, as in real 
ricepile experiments. This spirit is original from the philosophy of Ref. [7], although the 
model there is much more complicated. 

Even with the simplicity of its definition, Eqs. (1) and (2), the Oslo model gives rise 
to an astonishing complex behavior. As it is usual in SOC systems, it shows a power-law 
distribution of avalanche sizes [5,8] and avalanche durations [8], signaling the existence of no 
characteristic scales for the avalanche process. But also it has been shown that the transport 
of the grains through the pile is anomalous, in the sense that there is no normal diffusion 
but a power-law distributed transit time [5], spanning many orders of magnitude (just as it 
happens in the experiments, as we have mentioned). This has been explained by the fact 
that the time that a grain is trapped at a given position is also broadly distributed [9], 
which has in its turn being related to a kind of skewed fractional Brownian motion for the 
variations of the height, in the antipersistent case [10]. Moreover, the distances traveled by 
the grains during an avalanche turn out to be Levy flights [9], i.e., again scale free, despite 
the nearest-neighbor rules. Additionally, the time fluctuations of the profile scale with a 
roughness exponent that is in good concordance with the experiments [5,6]. In fact, the 
exponents of all previous magnitudes can be related by several scaling laws. 

Further, the time sequence of the transit times shows a clear multifractal spectrum [11], 
whereas the time sequence of the mean slope displays 1/ / noise [12], in contrast to the BTW 
model [2]. The model also allows one to study the transition from intermittent behavior 
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to continuous flow, just increasing the driving rate and breaking the time-scale separation 
[13]. Very recently it has been shown that a small damage performed in the system does not 
spread, and therefore the sensitivity of the model to the initial conditions is quite different 
to that of chaotic systems and to what was expected from systems "at the edge of chaos" 
[14]. On the other hand, there exists an exact mapping between this model an a model 
of interface depinning, which establishes the existence of a wide universality class for these 
non-equilibrium systems related to the quenched Edwards- Wilkinson equation [15-18]. In 
a previous paper, we have also signaled the similarities between the Oslo model and the 
recurrence of real earthquakes [19]. For all these reasons we can consider the Oslo model as 
the analog of the Ising model for slowly driven complex systems. 

In spite of these fascinating properties, our understanding of the Oslo model comes mainly 
from computer simulations and some scaling arguments; no analytical solution exists nor 
seems possible in the near future (as it is the general case for nonequilibrium systems). 
Hence, the exact enumeration of the number of states in the attractor for this model, per- 
formed by Chua and Christensen, is very remarkable [20]. They found that this number 
increases exponentially with system size as 

Ar aG L + a^G- L 

N A = ^ , (4) 

with G the golden mean, G = (3 + V5)/2 ~ 2.6, and a = (1 + y/E)/2 ~ 1.6. 

In this article, we are going to derive exact expressions for the transition probabilities 
between states in the Oslo model; using the results for a system of size L-l we will get these 
probabilities in a system of size L. The increase in just one unit of the size of the system 
leads to an enormous increase in the complexity of the resulting equations; we have a kind 
of machine for building complexity. With the transition probabilities it is possible to obtain 
the stationary occupation distribution, which is the probability with which a state is visited 
in the asymptotic regime. We will get that, in contrast to other SOC models, the states 
in the attractor are not equally likely; rather, the range of occupation probabilities varies 
dramatically in many orders of magnitude. Once the stationary occupation distribution is 



known, several other quantities, directly related with the profile of the pile, as the mean-slope 
distribution and the avalanche-size distribution can be obtained. (We have very recently 
become aware that D. Dhar has undertaken an analysis of precisely the same problem [21]; 
nevertheless, his approach is entirely different to ours and both works can be considered as 
complementary of each other.) 



II. SOME PROPERTIES OF THE MODEL 

Two types of states, or configurations, are possible in the system: unstable states, with 
at least one local slope value above its local threshold, z(x) > z th (x), and metastable states, 
where all the slopes are below threshold, z(x) < z t h(x) Vx. Unstable states evolve by means 
of avalanches towards metastable states, but the addition of a new grain can make metastable 
configurations become unstable again, and so on. 

We are only interested on metastable states. If we assume that the initial slopes z(x) 
are not negative [22], the possible stable values for this variable are 0,1, and 2 (values > 3 
are always unstable); therefore, in a system of size L a metastable configuration will be fully 
specified by an L-vector whose components are the slope values 0, 1, or 2. For instance, we 
can consider a state s to be s = {210 ... 1} which means z(l) = 2, z(2) = 1, z(3) = 0, and 
so on up to z(L) = 1 (alternatively, the metastable state s can be viewed as an integer with 
L digits expressed in base 3). Notice that we do not need to keep track of the thresholds 
z t h(x)] this is so because the dynamics can be described in an alternative but equivalent 
way, using the following rule: 



if z(x) has just changed its value to: < 



1, or less, no toppling, 

2, toppling with probability p, (5) 

3, or greater, toppling. 

We stress that this rule has to be applied only when the site at x has received one (or more) 
units of slope or has toppled in the previous avalanche (fast) time step. (This works because 
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we have only denned two possible values for the threshold, with three values the situation 
would be more delicate.) 

A very useful property in order to study the evolution of the system will be the Abelian 
symmetry, first considered in sandpile models by Dhar [23]. It states that the order in 
which units of slope are added and sites over threshold topple does not matter for the 
final configuration of the pile; therefore, we will be allowed to topple the sites in the most 
convenient sequence to keep the process manageable in the calculations. The demonstration 
of this property in our case is similar to that in Refs. [23-25] but taking into account that 
we have evolving thresholds. If we consider two sites x and y that are unstable it is easy to 
see that we get the same state no matter which one topples first, since after the toppling 
of x, site y will still remain unstable, and the quantity any toppling site is reset and the 
quantity transferred to the neighbors will be the same, independently of the order. The same 
reasoning can be extended to more than two over-threshold sites. But this is so only if the 
random thresholds for sites x and y are equally chosen in each possible sequence of topplings, 
that is, we need a predefined sequence of thresholds at each site, or, from a computational 
point of view, instead of having a single random number generator, a different generator 
must be used for each site. From a similar reasoning as before, the addition of grains (or 
slope) at x — 1 commutes with the toppling of any unstable site. 

On the other hand, the evolution of the pile can be described by means of a finite Markov 
chain; indeed, the probability that a given state transforms into another state depends only 
on the two states, and not on the previous history, thanks to rule (4). In particular, the 
probability that a metastable state i evolves to a new metastable state j after the addition 
of one grain and at the end of the corresponding avalanche is independent on the previous 
states of the pile, and can be obtained by means of the unstable states that separate the 
states i and j. These probabilities constitute a matrix that will be referred to as W, with 
elements Wij. Probability theory imposes that Wy > and that the files of the matrix are 
normalized to one (in the probabilistic sense), i.e., J2j = 1. A matrix with this properties 
constitutes a stochastic matrix. (Further, as W is constant through the time evolution, we 
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are dealing with a homogeneous Markov chain.) 

III. CHARACTERIZATION OF THE ATTRACTOR AND EXISTENCE OF A 
UNIQUE STATIONARY DISTRIBUTION 

Some results from the theory of Markov chains can be applied at this point. To get 
the stationary properties of the pile it will be crucial to have a well defined stationary 
occupation distribution; this quantity gives the probability with which every state is visited 
in the asymptotic regime, that is, in the attractor, and it is represented by a vector where 
each component corresponds to a configuration of the system. The stationary occupation 
distribution is simply referred to as the stationary distribution in Markov-chain theory, but 
here we are interested in many other probability distributions in the stationary for 
instance that for the avalanche sizes. (Another common name is ergodic distribution.) 

For completeness, let us explain that the attractor can be defined as the set of recurrent, 
or persistent states, being these the states for which the return probability is exactly one. 
More precisely, if a state is visited at some time, there is a probability one that it will be 
visited again in the future. In contrast, for transient states this probability is smaller than 
one, or even zero. 

At this point it is convenient to use graph theory to represent a Markov chain: the 
transition probability matrix W defines a graph where nodes % and j are directly connected 
if Wij 7^ 0, otherwise, there is no direct link between i and j; that is, we have the graph of 
the possible transitions in one (slow) time step. 

The existence of a unique stationary occupation distribution, independent on the initial 
conditions, is guaranteed if the graph associated to the matrix W has only one non-periodic 
final class [26], (this is also a necessary condition). A final class is a strongly connected 
component whose elements have no transitions to elements outside the class. A strongly 
connected component is a part of the graph in which any pair of nodes, or states, can be 
connected in both directions; in other words, i and / stay in a circuit and it is possible to 
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reach state / starting from % and vice versa. A final class represents then nothing else than 
an attractor in which the system can settle after a transient period. The periodicity of a 
strongly connected component is the greatest common divisor of the length of their circuits; 
if this number is one the graph is non-periodic. This is for instance the case in the presence 
of loops (circuits with just one element, that is, Wa ^ 0, for some i). 

Let us see which states of the pile constitute the attractor, or final class. We have found 
simpler to consider the connections between two states by means of the steepest metastable 
state, which is {22 ... 2} (i.e., the one with z(x) = 2, Vx), and then show which states lead 
to the steepest state and which ones result from it. 

In fact, all states can lead to the steepest state. To show this, we add grains and let the 
sites toppling depending on their local thresholds, but after every toppling we assume that 
the maximum threshold z th {x) = 2 is always selected; then, we are essentially in the same 
situation as in the one-dimensional BTW model: every column in the pile grows to reach 
the steepest profile. In this way, adding enough grains we will end in the steepest state. 
This is more easily seen applying the Abelian symmetry: we first add grains to built the 
first column, until it reaches the desired height, h(l) = 2L, then we add more grains and let 
them topple to the second column until it reaches h(2) = 2(L — 1), and so on. 

This is enough to ensure that there is only one final class, although the only thing we 
now about it is that it includes the steepest state. If we continue with the characterization 
of the final class we will simply get the attractor studied in Ref. [20]. 

On the contrary to the previous situation, not every state is reachable from the steepest 
state. Consider first final states without zeros, i.e., Zf(x) = 1 or 2, only, Vx. One way to get 
these states is the following: after the addition of the first grain, which crosses the whole 
pile arriving to the exit, we fix the thresholds to the slopes of the desired final configuration, 
z th( x ) = z f( x ), Vrr. Then, we apply the Abelian property and start the toppling process 
of the remaining grains from the rightmost site x = L, emptying out this column [from 
z(L) = 2 to z(L) = Zf(L)]; after this, we take the next column to the left, L — 1, and let 
every extra grain topple until it leaves the pile. Repeating the same procedure we end in 
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the desired state (which is therefore reachable after just one avalanche). Basically, we are in 
a one-dimensional BTW-like situation again, where the pile tends to a state z(x) = z th (x). 

When there are zeros in the configuration, z(x) = for some x, we cannot apply this 
trick since thresholds are defined as larger than zero. In fact, Chua and Christensen [20] 
have noticed that these states do not necessarily belong to the attractor: they argue that 
zeros have to be compensated by twos (i.e., sites with slope z(x) = 2); to be precise, they 
show that a necessary condition to belong to the attractor is that for each zero-slope site in 
the configuration there must be at least one two-slope site to the right, before the next zero 
or before the exit. 

Let us see that Chua and Christensen's condition is also sufficient to belong to the 
attractor: starting from the steepest state, the first zero appears when (after a number of 
topplings) a site with z(x) = 2 topples (if z t h(x) = 1) and the following site has z(x + l) = 1. 
Application of the rules gives then z(x) = and z(x + 1) = 2; if z th (x + 1) = 2 this site 
does not topple and the configuration can be metastable. Now that a zero exists, the same 
process can be repeated but with the zero-slope site receiving a grain; that is, we can have 
z(x — 1) = 2 and z(x) = 0, if z th (x — 1) = 1 we get z(x — 1) = 0, z(x) = 1, and z(x + 1) = 2; 
this means that the zero can move to the left, but is somehow associated to the existence 
of a site with slope two, and this slope two cannot disappear if the zero exists. In order to 
topple, site x + 1 would need the addition of one grain, but grains come from the left and 
cannot reach x + 1 except if the zero is removed, i.e., any grain coming from the left would 
encounter the zero slope and by the rules of the model would stick there. 

In general, to get a configuration {. . . 01 . . . 12 . . .} from the steepest state we go first to 
the configuration {...11. ..11...} (the same but replacing the and 2 by two l's) with the 
thresholds equaling the slopes. So, any new added grain will travel the whole pile up to the 
exit. After this, for the site that has to have slope 2 we make its threshold 2; an additional 
incoming grain will make the slope at this site equal to two and that of the preceding site 
equal to zero; successive incoming grains will move the position of the zero to the desired 
site. When there are more zero-two pairs in the final configuration the generalization is 
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straightforward if we start the previous procedure from the right. 

This demonstrates that the condition that each zero must have a two to the right is a 
sufficient condition to belong to the attractor. But further, the previous reasoning shows 
that states violating this condition are not accessible from the steepest state, nor from any 
other state which verifies the condition. So, the condition is necessary and sufficient, and the 
attractor proposed by Chua and Christensen constitutes the only final class of the system. 

Finally, it is easy to show the existence of loops in the final class: any state without zeros 
can return to the same state after the addition of one grain if this grain travels through the 
whole pile and does not induce any other grain to topple. For this we need that sites with 
slope one have also thresholds equal to one and sites with slope two keep their thresholds 
equal to two after toppling. (The probability of this is p n g L -™ ; where n is the number 
of sites with z(x) = 1.) This ensures the non-periodicity of the graph and completes the 
demonstration of the existence of a unique stationary distribution of state occupation. In a 
case like this, the Markov chain is said to be regular. 

Now that we know that there exists a stationary distribution, how do we obtain it? 
Notice that all that we have already learn about the system has been accomplished without 
explicit knowledge of the transition probabilities W; the relevant issue was if the transition 
was possible, Wij ^ 0, or not. However, to calculate the stationary distribution and proceed 
further we need the calculation of the matrix elements. 

IV. CALCULATION OF THE TRANSITION PROBABILITIES 

It is possible to derive the transition probabilities between the metastable states in a pile 
of size L as a function of the transition probabilities for a size L — 1. Since we have fully 
characterized the attractor [20] , we restrict the calculation only to these states, for the sake 
of conciseness. Note that although the number of states in the attractor [Eq. (4)] becomes 
astronomically large for the usual values of L in the simulations, it constitutes a drastic 
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reduction in front of the number of metastable states, which is 3 L ; i.e., 



/ \ 

V Q 



N A ~ 2.6 L < 3 L , 

for L large. 

Starting from L = 1 there are only two possible states in the attractor, {1} and {2} (state 
{0} is clearly transient since the coordinate x — 1 corresponds already to the boundary, 
where the toppling rules for the slope are special). We will label these two states as 1 and 
2 (in this case the label is straightforward, but not for larger L). Applying the rules of the 
model we easily get, 

where the superindex (1) stresses that we are dealing with a system of size L — 1. 

We can now consider L = 2 and look at the different states there, for instance, {11}. 
What happens when we add a grain at x — 1? There is a probability p that the origin 
topples, if not, we end in a state {21} with a probability q. If the origin topples, the grain 
jumps to position x = 2 and there the problem reduces to an L = 1 problem, for which we 
know the transition probabilities. In this case, we have to apply the transitions of state 1 
(which remember are defined taking into account that one grain is added to this state at its 
leftmost position, x = 2 now); as these transition probabilities are Wyy and W$ , we have 
for the state {11} the probabilities: 

q to go to state {z(l)z(2)} = {21}, 

vW$ = p 2 to return to {11}, 

J (6) 

p W[l ] = pq to go to {02}, 
for any other final state. 

This simple example shows how to reduce the problem from L to L — 1. In general, we 
will refer to the L system as the pile, or just the system, and the L — 1 pile, defined by 
x — 2,3 ... L, will be the subsystem or subpile. In the same way we will talk about states 
of the system and about substates when referring to the subsystem. 



W 



(2) 

{ll}{*(l)z(2)} 
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Although the previous case illustrates the basic idea, there appears an extra complication 
if the height at the origin, h(l), is larger, which is that the origin can topple more than once 
if the avalanche in the L — l subpile is big enough to leave the origin with a too large slope. 
In this case one can apply the Abelian property: let us topple first the subpile (just using 
the transition probabilities W^ L ~^ that we know) and then, at the end, let the origin topple. 
Of course, this gives rise to an iterative process, where the procedure has to be applied as 
many times as the origin topples, which is hi(l) — hp{l) + 1, / and F referring to the initial 
and final states. The sequence is: first the origin topples, then, the L — l subpile to reach 
equilibrium, then the origin again (if needed), then the subpile, and so on. 

Applying the previous argument to a general system of size L one can find the rules for 
the transition probabilities. It is convenient to define a variable Q as 

L 

Q = h(l) -L = 5>(x) -L- 

x=l 

we have Q — 0, 1 . . .L in the attractor [20]. The equations for the elements of will 
depend on AQ = Qi — Qf, that is, the difference of Q between the initial state and the final 
state (defined here in the opposite way as usual); this is so because the number of topplings 
at x = 1 is AQ + 1. The rules are given below and refer to the transition probabilities 
between an initial state / with a value of z(l) = zi and L — l subpile state i and a final 
state F with z(l) = zp and substate /, that is, a transition (zi,i) — * (zF,f). Note how 
an L— state is completely characterized by the value of z(l) and the state (substate) of the 
L — l subsystem. As with every toppling of the origin Q decreases in one unit, we will use 
the following relation to calculate z(\), 

z{l) =Q-Q' + 1, 

where Q' refers to the Q- value of the L — l subsystem. In general, Q' s will denote the Q- value 
of substate s in the L — l subsystem. T(z) will give the probability that a site topples for a 
given value of z [0, p, 1 for z < 1, 2, > 3 respectively, see Eq. (5)]. Therefore, the argument 
of T is the value of z(l) calculated after a number of topplings. With all these definitions 
the rules turn out to be, 
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w, 



(L) 

(zi,i)(z F ,f) 



= < 



[1 _ T { ZI + l)]S if 

T( Zl + l)^" 1 ^! - T(Qj - Q' f + 1)] 



if AQ < -2, 
if AQ = -1, 
if AQ = 0, 



T( ZI + 1) Ej W^TiQj - Q) + - T(Qj - Q' f )} if AQ = 1, 

+ 1) wg-VTiQj - Q' 3 + 1) 



if AQ = 2. 

(7) 



The case AQ < —2 is impossible since we would need adding more than one grain at x — 1 
(and that some of them would not topple) to reach the corresponding value of Q. For 
AQ = — 1 the only possibility is that the origin does not topple, then zj increases in one 
unit and has to be stable, and the substate does not change. For the rest of cases we can 
write 



W, 



(L) 



(z I ,i)(z F J) 



T( Zl + l) Ej w t f- 1] T(Qi - Q] + i) Ek W$- 1] T(Qi - Q'u) 



(8) 



Ez W% ~ 1] T(Qj - Q\ - 1) • • • Wi^T(Qj -Q' v -AQ + 2) 
wi L f - l) \l - T(Qj - Q' f - (AQ - 1))] if AQ > 0. 

The general idea is that there is a probability T(zj + 1) that the origin topples after the 
addition of one grain, then, there is a probability W^~^ to go to a substate j; with this 
substate there is a probability that the origin topples given by T again, and then, one goes 
from j to k, from here to I, etc., following all the possible paths that end in /, whose origin 
has a probability 1 — T(z F ) to be stable, with z F = Qp — Q'f + 1 = Qi — Q'f — AQ + 1. 

The previous equations for the elements of look like a matrix product but with 

different matrices and a different number of factors for different components. Nevertheless 
in matrix form they can be written as 



W, 



(zi,i)(z F ,f) 



T( Zl + l)[l-T(z F )} 



1=Qi 



if AQ > 0, 



if 



where the index Q is assumed to decrease in the product, [ ],/ denotes that we take the 
element of the matrix in the % file and / row and T is a diagonal matrix whose elements can 
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only be l,p, or 0, more precisely, 



[T Q ]. k = T(Q - QJ + l)5 jk . 

We stress that the rules are valid for any pair of metastable states, although we will con- 
centrate on states in the attractor. 
These equations for L = 2 yield 
/ 



W (2) 












1 





pq 


p 2 


q 








2 

p q 


p3 


pq 


Q 





p 3 q 


4 

P 


2 

p q 


pq 


q 



(9) 



/ 



\p 3 q(l + q) p\l + q) p 2 q(l+q) pq(l + q) q 2 
where states are ordered as {02}, {11}, {21}, {12}, {22}. Iterating the rules it is possible, 
although laborious, to generate the matrices for successive L. Although the matrix for 
L = 2 looks simple, the corresponding matrices as L increases are getting more and more 
complicated. 

Figure 1 shows, for L = 7 and 8, and p — 1/2, the probability density H(K) of the 
number of nonzero elements for each file of W, that is, the distribution of the number of 
states directly accessible for a state in the attractor, or, in the language of networks, the 
out-degree distribution of the phase space. There seems to be two kinds of states, one group 
has few connections and the other one a large number of them; nevertheless the system 
size is too small to be conclusive. In contrast, the in-degree distribution (not shown in the 
plot) looks rather uniform. (We will always use the letter H to denote probability densities, 
although it will correspond to different functions depending on the argument.) 

In Fig. 2 we illustrate the enormous variation in the values of the transition probabilities: 
the probability density H(W) that Wij takes a certain value W is shown for L = 7 and 
p = 1/2, spanning about 16 orders of magnitude. A power law with exponent minus one 
approximates well this behavior. Curiously, a similar result has been found in networks 
describing correlations between earthquakes [27]. 
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V. OBTAINING OF THE STATIONARY OCCUPATION DISTRIBUTION 



To get the evolution of the system one has to start with a distribution of initial states P , 
which remember has to be an Na— dimensional vector [see Eq. (4)], giving the occupation 
probability of any of the Na states in the attractor. This initial distribution can be a delta 
[for instance, for L = 2, starting always with {11}, i.e., P = (0,1,0,0,0)] or not. The 
distribution of states after the first (slow) time step, i.e., after the addition of just one grain, 
is obtained as P\ — Pq • W. To obtain the distribution of states in the next time step we 
have to multiply again by W and so on; the powers of W give therefore the evolution of the 
system. 

Let us call the vector representing the stationary distribution D, which of course is also 
a vector on Na dimensions, with each component D(s) giving the probability that after a 
long enough time the system is in state s. The evolution of the occupation probabilities in 
the attractor is also obtained multiplying the row vector D by the matrix W, but as D is the 
stationary distribution, it must be invariant under such operation, i.e., D = D • W, which 
simply means that D is a left eigenvector of W with eigenvalue equal to one. Regularity 
(which was demonstrated in Sec. Ill) ensures that this eigenvector is unique. 

Moreover, a direct consequence of regularity, which is also a sufficient condition for it, is 
that 



M = lim W l = 



' D ' 



D 



(10) 



which means that the transition matrix corresponding to n steps has all its files equal to 
the stationary distribution D, asymptotically. Indeed, it is trivial to show that for any 
distribution P, with J2i P(i) = 1, we have P ■ M. = D. On the opposite side, if any P tends 
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to D we can take P = (1, ... 0) to get the first file of matrix M (by multiplication), which 
must be equal to D, and the same can be done for any other file of M. 

If we consider the case L = 2 we realize that M. = W 3 , that is, we obtain a matrix 
whose files are all the same; this implies that the asymptotic result is reached in just three 
time steps, since M = W ■ M. = W 2 ■ M, etc. But further, it turns out that the stationary 
distribution D, given by the files of Ai, is the last file of W, if the states are ordered by 
increasing Q, so, 

D = ( ptq(l + q), p\l + q), p 2 q(l + q), pq(l + q), q 2 ) ; 

that is, the transition probabilities of the steepest state {22} give the occupation probabilities 
in the stationary regime. In other words, the unique eigenvector of W with eigenvalue equal 
to one is just its last file. 

We have verified that this result is general for larger L, although the necessary number 
of powers increases with L (note that for L = 1 this is already accomplished at the first time 
step). Dhar [21] has beautifully demonstrated this result using the properties of an operator 
algebra. The idea behind this is simple: we can realize that it is equivalent to add L(L + 1) 
grains to the flattest state {00 ... 0} (which is not in the attractor) than to add just one 
grain to the steepest state {22 ... 2}. Let us see why. The number of grains which separate 
both profiles is precisely L(L + 1), so, by application of the Abelian property, we add this 
number of grains to {00 ... 0} and let them topple to reach the profile corresponding to the 
steepest state (the probability of this toppling process is exactly one); after reaching this 
state we can continue the toppling process, but we are already in the same situation that 
results from adding just one grain to the steepest state. So we get the same configurations 
with the same probabilities in both cases. In fact this result is not only true for the flattest 
state, but for any other state, the only difference is that we will have some extra grains: no 
problem, they topple until they leave the pile. Therefore we can write 

(0...0,1)-W = P -W L(L+1) , 
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VP = (0 . . . 0, 1, . . . 0). (This expresses something that we already proved in Sec. Ill, 
which is that the steepest state is reachable from any configuration, just adding enough 
grains.) Note now that the addition of one extra grain changes nothing in each case, this 
grain will topple until the exit, so 



which means that, indeed, the transitions from the steepest state coincide with the stationary 
occupation distribution. Dhar has also noted that a more restrictive condition holds for the 
states in the attractor. The state there with less grains is {11 ... 1}; the difference in number 
of grains between the steepest state and this one is L(L+l)/2, which can replace the previous 
value L(L + 1), for any state in the attractor. 

In Fig. 2 we also include the probability density H(D) that -D(s) takes a given value for 
L = 7 and p = 1/2, showing a behavior very similar to the density of transition probabilities; 
this is a broad distribution across 13 orders of magnitude, close to a power law with exponent 
minus one. This means that the occupation of the phase space (i.e., the space of all possible 
configurations) is enormously heterogeneous, at variance with the BTW model [23]. 

Figure 3 shows D(s) as a function of s, where the states s are ordered in terms of 
decreasing D(s). In fact, the form of D(s) in this plot is related to H(D), just by identifying 
s/N A as the probability that the occupation probability is larger than (or equal to) a certain 
value D(s), that is, as the survivor function of the random variable D(s). Therefore, the 
density H(D) will be (as usual) the derivative of this survivor function, multiplied by — 1, 
or 



P ■ W L(L+1) = P ■ W L{L+1)+1 . 



As this equality holds for any vector of the basis, we can write 



yyi(L+l) = yyL(L+l)+l = ^ 



and so, from the first equation we get 



(0 . . . 0, 1) • W = P ■ M = D, VP , 
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A power law with minus one exponent for H(D) yields an exponentially decreasing D(s), in 
agreement with the plot. 



VI. CALCULATION OF THE DISTRIBUTIONS OF MEAN SLOPES AND 

AVALANCHE SIZES 

From the values of the stationary distribution of the occupation of the states and their 
transition probabilities, D(s) and W, it is possible to calculate many things in the stationary 
state, (i.e., in the attractor). The first one is the (stationary) distribution of Q, f(Q), 



for Q = 0,1,2. 

The distribution f(Q) is in fact the distribution of heights at the origin, f(h(l)), and 
it is also directly related to the distribution of mean slopes, f(z), with z = J2x=i z ( x )/L = 
h(l)/L, since Q = h(l) — L = L(z — 1). Consequently, f(Q) defines the active zone width, 
which can be obtained as the standard deviation of this distribution. (For simplicity, we 
have used the same symbol / for all the distributions, although obviously they are not the 
same function.) 

The distribution for the value of (Q - (Q))/L* = (h(l) - (/i(l)))/L* = (z - (z))L l ~* is 
displayed in Fig. 4 for several values of L. For x, we take the value proposed in Ref. [13], 
X — 0.24. Note how all the discrete distributions collapse onto a single continuous curve 
under rescaling, which is close to Gaussian, though slightly skewed. These exact results for 
small L are in total agreement with the findings of computer simulations. 

The avalanche size distribution f(S) (in the attractor) is not difficult to calculate knowing 
D(s) and W. We can write 



HQ) 



E D(*)=T, D (*)*Q.Q- 



Vs s.t. Q S =Q Vs 



For example, for L = 2 we get 




f(S) = Y,P(S/s)D(s), 
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where S is the avalanche size, s the state of the pile, and p(S/s) is the conditional probability 
of having an avalanche of size S starting from a state s. For this term we have, 

P (S/s)= £ W aj = ^W aj 5 Saj 8, 

Vj s.t. S sj =S Vj 

where S S j is the size of the avalanche triggered in the transition from s to j and can be 
calculated as 

L 

s sj = J2( h s( x ) - h i( x ))( L -x + l) + L, 

x=l 

which is essentially the profile difference times the distance to the exit, plus the contribution 
of the added grain. Therefore we have, 

Vs Vj 

The corresponding distribution calculated in this manner for L = 8 and p — 1/2 appears 
in Fig. 5, where it is compared with the result obtained from computer simulations. Note 
how even for such a small system there are avalanches with probability smaller than 10~ 22 . 

VII. DISCUSSION 

Just to summarize, we have shown that the conditions about the states put forward in 
Ref. [20] are necessary and sufficient to define the attractor, which moreover is unique. In 
addition, its nonperiodical character allows the existence of a single occupation distribution 
in the stationary limit. The Abelian property enables the calculation of the transition 
probabilities between the different states, just from the decomposition of a system of size 
L into its leftmost column and an L — 1 subsystem and by using an iterative toppling 
procedure; this allows to explore the network of connections in phase space. Unexpectedly, 
the stationary occupation distribution turns out to be the last file of the transition matrix, 
i.e., that corresponding to the transitions of the steepest state. Both the transitions between 
states and their occupations can take values in an extremely large range of probabilities, 
setting a clear difference with other SOC systems. These calculations are exact for the 
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system sizes involved, and could be performed for any L, in principle. In practice, we 
have strong limitations, as it is explained below. Finally, with these quantities we can also 
derive the form of the fluctuations of the profile and the avalanche size distribution, for the 
corresponding value of L. 

In fact, the knowledge of the transition probability matrix and the stationary distribution 
allows the calculation of any property related to the profile of the pile, as the ones we have 
just mentioned or the dissipated-energy distribution. But there are other properties that 
do not only depend on the profile but also on the dynamics which leads from one profile 
to another, for instance, the avalanche duration: knowing the profiles is not enough to 
calculate this quantity; in fact, with the same initial and final states the duration is not 
uniquely defined, i.e., there are many ways to end in the same state with a different number 
of avalanche time steps, depending on the actual sequence of topplings. So, despite the 
usefulness of the Abelian property, it does not allow to calculate dynamical magnitudes. 
This could be solved in principle with a parallel calculation stating the probability that a 
given transition would involve a given time; nevertheless, it would be complicated. 

Also, the profiles, that is, the configurations defined in terms of slopes (or heights), do 
not retain any information regarding individual grains. Grains are essentially treated as 
undistinguishable whereas for the calculation of transit times or flight lengths one needs 
distinguishable particles. This is another limitation of the current method difficult to over- 
come. 

On the other hand, the equations obtained for the transition probabilities have the 
advantage that one does not have to simulate the system, and therefore one obtains exact 
probability distributions. However, the enormous number of states in the attractor (which 
increases as 2.6 L ) makes impossible any calculation by hand beyond the smallest values of 
L. Symbolic computer calculations, performed by MAPLE or similar programs, have also 
severe limitations of size, and even with numeric computations the system size is limited to 
about L < 10 (for L = 10 there are more than 10 4 recurrent states, which would require 
10 8 matrix elements). Also, the number of different matrix products is large. Although 
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the results we obtain for small systems are interesting enough and representative of the 
complexity of the system, it would be nice to have access to supercomputers to increase the 
capabilities to manage larger system sizes. 

To conclude, it is interesting to point out that having exact analytical results for some 
problem is not a synonymous of understanding; in this case we can obtain exact expressions 
for the transition probabilities, the stationary distribution, or the avalanche size distribution, 
but from these formulas still it is not clear which are the properties of the system for large L. 
For instance, for the avalanche size distribution we can get very complicated exact equations, 
but we cannot show that they tend to a power law in the asymptotic case. Nevertheless, in 
our opinion, the results in this paper imply a remarkable progress in our picture of complex 
systems. 
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FIGURES 
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FIG. 1. Probability density that the number of states directly accessible in phase space (or 
out-degree distribution) for a state in the attractor takes a value equal to K, for L = 7 and 8 and 
p = 1/2. For K not very large, H{K) could be a power law, but larger values of L are needed to 
be more sure. The histogram is calculated with exponentially increasing bins. 
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FIG. 2. Probability densities H(W) and H(D) that some transition probability takes a value 
W, and that some state s has stationary occupation probability D(s) = D, for L = 7 and p = 1/2. 



25 




200 400 600 800 1000 1200 1400 1600 

s (by decreasing D(s)) 



FIG. 3. Stationary probability of occupation D(s) for each state for L = 7 and 8, and p = 1/2. 
The states are ordered in decreasing probability. Straight lines would indicate an exponential decay 
with s. 
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FIG. 4. Stationary distribution of Q (or of heights at the top, h(l)), centered by the mean and 
scaled by L 0,24 . System size ranges from L = 1 to L = 8, and p = 1/2. Notice how all the discrete 
distributions for each system size conspire to give a smooth curve, which is close to Gaussian, 
although slightly skewed. 
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FIG. 5. Stationary distribution of avalanche sizes for L = 8 and p = 1/2 from our exact 
procedure and from simulations of the Oslo model. 
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